Measuring emission coordinates in a pulsar-based relativistic positioning system 
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A relativistic deep space positioning system has been proposed using four or more pulsars with 
stable repetition rates. (Each pulsar emits pulses at a fixed repetition period in its rest frame.) 
The positioning system uses the fact that an event in spacetime can be fully described by emission 
coordinates: the proper emission time of each pulse measured at the event. The proper emission time 
of each pulse from four different pulsars — interpolated as necessary — provides the four spacetime 
coordinates of the reception event in the emission coordinate system. If more than four pulsars are 
available, the redundancy can improve the accuracy of the determination and/or resolve degeneracies 
resulting from special geometrical arrangements of the sources and the event. 

We introduce a robust numerical approach to measure the emission coordinates of an event in 
any arbitrary spacetime geometry. Our approach uses a continuous solution of the eikonal equation 
describing the backward null cone from the event. The pulsar proper time at the instant the null 
cone intersects the pulsar world line is one of the four required coordinates. The process is complete 
(modulo degeneracies) when four pulsar world lines have been crossed by the light cone. 

The numerical method is applied in two different examples: measuring emission coordinates of 
an event in Minkowski spacetime using pulses from four pulsars stationary in the spacetime; and 
measuring emission coordinates of an event in Schwarzschild spacetime using pulses from four pulsars 
freely falling toward a static black hole. 

These numerical simulations are merely exploratory, but with improved resolution and compu- 
tational resources the method can be applied to more pertinent problems. For instance one could 
measure the emission coordinates, and therefore the trajectory, of the Earth. 

PACS numbers: 04.25.D- 95.10.Jk 



I. INTRODUCTION 

Pulsars — spinning neutron stars that emit directional 
electromagnetic radiation with an intriguingly stable 
period — in principle can be used as reliable interstellar 
lighthouses for spacecraft navigation in the Solar Sys- 
tem and beyond [IH^- Pulsar spacecraft navigation has 
begun to be experimentally investigated by observing X- 
ray pulsars ^3; from satellites, comparing these observa- 
tions to satellite ephemerides. Rovelli [5] suggested a 
method to construct this fully-relativistic universal coor- 
dinate system based on the proper emission times of the 
emissions from sources that he called "satellites". The 
concept of coordinates based on the emission times of 
pulses is identical in concept to the global positioning 
system (GPS). However, we treat a fully relativistic for- 
mulation, while the current GPS system treats relativis- 
tic effects only as perturbations from a Newtonian frame- 
work. A number of authors have developed the idea of 
a completely relativistic satellite positioning system; see 
[SHU] • Delva et al. [13] gave a recent review of this idea, 
including an extensive reference list. 

Consider four pulsars, in motion in space, broadcast- 
ing pulses at a constant rate as measured in their proper 
tinres. The intersections between the world lines of these 
pulsars and the past light cone of a reception event D? 
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give the proper emission times of the pulses from each 
pulsar that will be recorded at the event 3J. (We elabo- 
rate on the definition of the pulsar proper time, and on 
the interpolation between pulses in Section |VI| below.) 
The event "R, then, can be described by the coordinates 
(t^ , t'^ , T^ , T*) called the emission coordinates, where r^ 
refers to the proper emission time of the pulse from pul- 
sar #1, T^ refers to the proper emission time of the pulse 
from pulsar #2, and so on. A spacecraft recording the 
proper emission times of pulses from the four pulsars will 
then be able to determine its coordinates — and therefore 
its trajectory — in spacetime. 

These emission coordinates can then in principle be 
converted into more conventional spacetime coordinates 
{t,x,y,z). This paper provides a proof of concept nu- 
merical demonstration of determining the emission coor- 
dinates, and converting these coordinates to more stan- 
dard spacetinre coordinates. 

Our method makes use of level set solutions of the 
eikonal equation describing the past light cone of the 
event 3?. This method was originally developed by 
Caveny, Anderson, and Matzner to track black hole event 
horizons in computational simulation of black hole inter- 
actions [TJj , and is similar to approaches in Refs. [TS] and 
|16) . In the present paper, we show that the same numer- 
ical approach can address the problem of interconverting 
spacetime coordinates and the respective emission coor- 
dinates. This approach is complete in the sense that a 
single method can be used to measure these emission co- 
ordinates even when the observer at 3? — our spacecraft — 
and the pulsars are moving in an arbitrary spacetime ge- 



ometry. 

The problem of converting between conventional and 
emission coordinates naturally arises as one begins to 
develop an intuition for emission coordinates. It was 
treated extensively by Delva and Olympio [IT] ; they have 
in mind that the source is a navigational satellite in the 
Schwarzschild spacetime representing the Earth's grav- 
itational field. Our eikonal solution for the backward 
null cone of the reception event Jl adds a new method to 
determine the emission coordinates of Jl, in addition to 
those of dl]. 

Our method is given in a proof-of-principlc form, with 
moderate computational accuracy; we discuss means to 
improve its accuracy (Section VII below) . Its advantages 
are that it is general and robust in any spacetime (we 
give flat space and Schwarzschild examples); it involves 
no "shooting" or other iterative methods; it involves no 
approximations except discretization for computational 
integration; it builds up the entire past null cone; if the 
past null cone intersects the source world line, the in- 
tersection, and hence the emission coordinate, will be 
found. To emphasize to the generality of the method: the 
Minkowski and the Schwarzschild examples differ only in 
the metric used. The spinning Kerr spacetime could be 
treated similarly by inserting the Kerr metric instead. 
And the method will straightforwardly work with a met- 
ric given only computationally, for instance the result of 
a simulation of the gravitational field of a binary pulsar. 

This work is organized as follows: Section |ll] outlines 
the theoretical framework. Section Hill describes the im- 
plementation of the numerical description of the eikonal 
equation. Our numerical approach uses a second-order 
artificial viscosity term, as in Caveny et al. |14j . Section 
|IV| presents the results of a numerical simulation where 
all pulsars are placed stationary in flat Minkowski space. 
Section|V]presents the results of a curved space numerical 
simulation where the pulsars are freely falling toward a 
Schwarzschild black hole. In both Section ITVl and Section 
[yj we actually measure the emission times from five pul- 
sars. This allows us to construct a configuration which 
is easy to plot: four pulsars in the conventional coordi- 
nate plane z = 0, a configuration which however has an 
obvious degeneracy between ±z when determining the 
emission coordinates (see Figures Il| and [2] below) . The 
fifth pulsar is chosen out of the z = plane, and emission 
coordinates based on pulses from the fifth and any three 
of the first four pulsars are free of the z — degener- 
acy. But such a system is difficult to plot graphically, 
and we do not try. The multiple coordinate systems 
constructed in this way are related to one another by 
flnite coordinate transformations. Section IVII discusses 
this and the continuous gauges in which the emission co- 
ordinates are measured, and some improvements that are 
required when constructing a more practical pulsar-based 



positioning system. Section VII discusses future improve- 
ments to the current code that will improve the code's 
numerical accuracy. Conclusions are presented in Section 

Eml 



We use geometrical units throughout, so the Newton 
constant G and the speed of light C are set equal to unity. 
We also use the Einstein summation convention so that 
repeated (one up, one down) indices are summed over 
their range. Greek lower case indices range and sum over 
{0, ..,3}; Latin lower-case indices range and sum over 
{1,2,3}. 



II. THE EIKONAL EQUATION 

See Ref. [13]. The world line of a photon can be de- 
scribed by the equation of motion. 



d /9L 

dr V9x" 



dL 
dx°' 



= 0, 



(1) 



where r is an affine parameter and x" = ^^. Since the 
Lagrangian of a null geodesic motion, L — ■kgapi'^i = 
has only kinetic terms, it is equal to the associated 
Hamiltonian — obtained using the Legendre transforma- 
tion: 



H = ^g^^PcPp = L. 



(2) 



We introduce the 3-1-1 Arnowitt-Deser-Misner vari- 
ables a, l3j and 7,^, 



7jj — Siji Pi 



9ti 



M'-9tt, 



(3) 



where indices on f3i are raised by 7*-' (the three- 
dimensional inverse of jij) and lowered by jij. Equation 
([3]) implies 
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(4) 



For a photon, which follows a null geodesic motion, 
the Hamiltonian has value zero and we can solve for pt 
to find 



Pt = 13' Pz ± a^-f^ipiPj. (5) 

The eikonal equation can be obtained by making sim- 
ple direct substitutions Pt ^ ^ = S^t and Pi — ^ ^ = 
S^t in Eq. pi (dropping the factor of 1/2): 

9"^S,^Sj - (6) 



which can be solved for Sj. Using Eqs. (IS]) and Q, we 
obtain the following symmetric hyperbolic partial differ- 
ential equation 



S,t^(3'S,,±aJr'S,,S,, = -H. 



(7) 



The bar is used to distinguish the Hamiltonian used here 
from the Hamiltonian in Eq. (l2|. iJ is homogeneous of 



degree 1 in S,i. The characteristic curves along which the 
level sets F of 5' are propagated are therefore 






(8) 



and 



= -|.(/^%±"v/^^^^)="S^ (9) 



dS,, d 



which are the null geodesic equations. The integral 
curves of the gradients of 5* and T are then also the null 
geodesies: 



da:'(A) 
dA 



g>a=5*"^,a = 5'''(A,x^(A)) 



(10) 



where A is the affine parameter in this case. 

Solutions to the eikonal equation fall into topologically 
equivalent classes. This means that any smooth function 
il]{S) topologically equivalent to S is also a solution. Note 
the following relations, 



*•<*)- IS -^<*)*- 



and 



^AS) 



(hl^dS_ 
dS~dt 



HS)s^t 



(11) 



(12) 



The above relations are true because Eq. (It]) is ho- 
mogeneous of degree 1 in momentum. The above results 
guarantee that a smoothly related initial data So ^ S'q = 
ip{So) have smoothly related solutions. 

The solutions of the eikonal equation also guarantee 
the equivalence of ingoing and outgoing solutions under 
time reversal. Referring to Eq. nh, propagation of data 
for S describing an ingoing or outgoing null surface is 
completely specified by: 

1. a definition of the direction of time, 

2. choices of a and /3*, and 

3. a choice of the sign of the root. 



III. NUMERICAL METHOD USING THE 
EIKONAL 

Our numerical method makes use of the time evolution 
equation for the solutions 5, 



S,t^-Hit,x\S,,), 



(13) 



which allows for relatively fast calculations while being 
sufficiently accurate. 

The eikonal equation, however, shows singular behav- 
ior and as described by Ehlers and Newman [T^, the 
eikonal equation generally breaks down on caustic and 



other sets. To address this problem, we make use of an 
explicit artificial viscosity term to control the appear- 
ances of such singularities. Adding the artificial viscosity 
at the level of the finite difference approximation corre- 
sponds to replacing the time evolution Eq. ( 13 1 with the 
evolution equation 



S^t^ey^S-H{t,x\S^, 



(14) 



where e — the artificial viscosity — is a small quantity of 
the order of h"^ {h denotes the resolution of the numeri- 
cal mesh) and V^ is any second-order, linear derivative 
operator. We use a second-order finite difference approx- 
imation to the Laplacian. 

The null surface F, at any given time level, can be ex- 
tracted from the level set section of the eikonal solution 
S, say 5 = 1. This problem of extraction is an inverse 
problem, since it requires that points {x,y,z) are found 
such that S{x,y,z) = 1. Nevertheless, a combination 
of ordinary bisection and interpolation method is suffi- 
cient to extract the approximate null surface F. In this 
method, the surface F can be represented in spherical co- 
ordinates {u{9 , (j)) , 9 , (j)) , where r = u{9,(j)) is the surface 
function for a given center c* contained within the surface 
f. 



IV. NUMERICAL SIMULATION IN 
MINKOWSKI SPACE 

The first application of our numerical method to mea- 
sure emission coordinates is done in the simplest con- 
figuration possible: the pulsars are stationary in flat 
Minkowski space and we are measuring the emission co- 
ordinates of an event 31 at the spatial origin (t, x, y, z) — 
(0.1,0,0,0). Although Minkowski space has no natural 
timescale, we may take the coordinates to have units of 
seconds. 

As described in Section II, defining the numerical sim- 
ulation requires the choices of the direction in time, a, 
/3* (determined from our Minkowski coordinates as a = 1 
and Pi = 0), and the sign of the root in Eq. (I?]). To mea- 
sure the emission coordinates — the intersections between 
the past light cone and the world lines of the pulsars — we 
need to choose the sign of the root to be negative so that 
the propagation of S describes an outgoing null surface 
when the direction of time is pointing to the past. 

The simulation is done in a three-dimensional compu- 
tational domain of N^ points with N = 361. The outer 
boundaries are located at [—2.5, +2.5] in the x,y^z di- 
rections. The resolution of this finite difference mesh is 
then h = 5/360 sec. A Courant-Friedrichs-Lewy factor 
of A = 1/4 with an iterated Crank Nicholson scheme [18] 
is used in the finite difference approximation of the time 
evolution equation ( 14 ) . Here the artificial viscosity pa- 
rameter e is set to be equal to h^ /16. We found that 
applying w 400 numerical evolution steps (stopping at 
the time slice t = —1.388 sec) was more than enough to 
measure the emission coordinates. 



In this simulation, the stationary pulsars are located (111 and (12). What Eq.(16) accomplishes is to smooth 



at: 

1. pulsar #1: (t, x, y, z) = (t, -0.50, 0, 0), 

2. pulsar #2: (t, x, y, z) = {t, 1.00, 0, 0), 

3. pulsar #3: (t, x, y, z) = (t, 0, -0.75, 0), 

4. pulsar #4: (t, x, y, z) = {t, 0, 1.25, 0), and 

5. pulsar #5: {t,x,y,z) = (t, 0.30, 0.40,0.50). 

The first four create a quadrilateral in the z = plane 
surrounding the observer at x = y — z — 0^ which will 
generally produce a ±z degeneracy when determining the 
emission coordinates. This configuration is easy to plot; 
see Figure [l] Including pulses from the fifth pulsar with 
those from three of the others create emission coordinates 
without the ±z degeneracy. Alternately one can combine 
data from more than four (e.g. five) pulsars by a kind 
of least squares fitting. That is the subject of work by 
Tarantola et al. [10] . We assume that the proper time of 
each pulsar is precisely the Minkowski coordinate time t; 
this choice is available only in Minkowski spacetime, but 
not in curved spacetime (see Sections |v] and VI I . 

In order to avoid the focal point singularity at the ver- 
tex D? of the null cone, we represent the event by positing 
a spherically symmetric null surface T centered at the 
origin with radius p = 0.1 sec at time t — sec. We thus 
assume (trivially correct in Minkowski space) that near 
the event under consideration, spacetime is sufficiently 
fiat so that the light cone is spherical around the event, 
and set the data a small amount of time (here 0.1 sec) in 
the past of the event we are coordinatizing. 

Data for the eikonal equation are set in the spherically 
symmetric form 



S(t = 0, x') = 1 + tanh 



(15) 



where Tc is the radius of the initial surface T and is equal 
to 0.1 sec in this case, while c controls the steepness of 
the hyperbolic tangent function. We set c to 0.1 sec in 
our experiment. 

After every 75 iterations, we reinitialize the data for 
the eikonal equation with a function that is similar to 
Eq. (pi). 



S'(t,a;') = 1 + tanh 



u{0, (j),i) — r 



(16) 



to ensure the smoothness of the data: reinitialization al- 
lows us to set viscosity parameter e to be arbitrarily small 
while avoiding the onset of singularity. Here c' denotes 
a new steepness of the function and our simulation uses 
c' — c. Recall that u{9, cj), t) is the surface function for a 
given center c' contained within the discrete null surface 

r. 

We know for a fact that the reinitialized solution S is 
also a solution to the eikonal equation by virtue of Eqs. 



the function S near the location of the null surface (where 
5 = 1) to produce data for a solution which are both 
analytic and smooth, and which describe the same null 
cone. 

For the resolution used in our simulations, we find no 
difference in behavior in the Minkowski case, whether or 
not this reinitialization is carried out. However, in the 
Schwarzschild curved spacetime case, the reinitialization 
is necessary, as discussed below. 

TableUshows the results of our measurements to deter- 
mine the emission coordinates of the event (t,x,y,z) — 
(0.1,0,0,0). 

Figure [T] depicts the intersections of the past light 
cone generated by our numerical simulation in the equa- 
tor {z = 0) with the world lines of pulsars #1 through 
#4, depicted by the four bold lines. The dots mark the 
events when each pulsar emits its pulse, which will later 
be recorded by the observer at the event 31 we are coor- 
dinatizing. 



Pulsar 


# 


r 


t 


X 


y 


z 


1 




-0.399 ±0.003 


-0.399 ±0.003 


-0.500 


0.000 


0.000 


2 




-0.903 ±0.003 


-0.903 ±0.003 


1.000 


0.000 


0.000 


3 




-0.653 ±0.003 


-0.653 ±0.003 


0.000 


-0.750 


0.000 


4 




-1.156 ±0.003 


-1.156 ±0.003 


0.000 


1.250 


0.000 


5 




-0.615 ±0.003 


-0.615 ±0.003 


0.300 


0.400 


0.500 



TABLE I. Results of emission coordinates in Minkowski 
spacetime of an event D? = (t,x,y,z) = (0.1,0,0,0). r de- 
notes the proper time of the pulsar when the pulsar world line 
intersects the observer's past light cone. Any four of the five 
proper times listed in the table constitute the emission coordi- 
nates in a particular emission coordinate system, of the event 
point 3?. The coordinates i, x, y, and z are the Minkowski 
space coordinates of the pulsars when the intersections occur. 

Comparing our results with analytical calculations, the 
errors in our results are mainly of the order of a few time 
resolutions of the finite difference mesh, dt ^ X ■ h ^ 
0.0035. Simple analytical calculations show that world 
lines of pulsars #1, #2, #3, #4 and #5 are expected 
to intersect the null surface at t — —0.400 sec, —0.900 
sec, —0.650 sec, —1.150 sec, and —0.607 sec, respectively. 
The errors can be attributed to the numerical accuracy 
of the simulations, mainly due to: (1) the interpolation 
routine that extracts the discrete surface F, and (2) the 
resolution of the underlying three dimensional grid. 



V. NUMERICAL SIMULATION IN 
SCHWARZSCHILD GEOMETRY 

Our numerical method has been successful in measur- 
ing the emission coordinates in Minkowski spacetime but 
in general spacetime is not flat. To investigate the change 
in the emission coordinates in a curved spacetime, we 



a(> 0) and /3i are: 




15 2 



.2 -1.5 



FIG. 1. Rays: plot of null surface f (the past light cone of 
K in Minkowski spacetime). Star marks the location of the 
event 3? = (t, x, y, z) — (0.1, 0, 0, 0) that we are coordinatizing. 
Bold lines: world lines of the four stationary z — pulsars in 
Minkowski spacetime. Dots mark the intersections between 
the null surface and the pulsars. The coordinates of the inter- 
sections are recorded in Table [I] Analytic methods are used 
to substitute for the null cone between K and the top of the 
computed (truncated) null cone. See the discussion in the 
text. 



evaluate the numerical method in a Schwarzschild geom- 
etry containing a stationary black hole of mass M = 0.25. 
We use the Eddington-Finkelstein [121 HHI coordinate 
system in describing the Schwarzschild geometry to avoid 
the coordinate singularity at areal distance r = 2M 
from the black hole. The standard Eddington-Finkelstein 
form of the Schwarzschild solution is centered at the 
spatial origin of the Eddington-Finkelstein coordinates. 
In order to maintain our event position at {x,y,z) = 
(0,0,0), we offset the black hole coordinates, putting it 
at (xo,2/o,zo) — (2-5,0,0) to provide a strong gravita- 
tional attraction at (0, 0, 0). The Schwarzschild metric in 
Eddington-Finkelstein coordinates (in Kerr-Schild form) 
is ini: 



2A'f, , 

gal3 = flajS H tQt/3 



(17) 



where r = ^J{x — xoY + (y ^ VoY + (-^^ ~ '^oY is the 
areal coordinate distance from the center of the black 
hole. Here riap = diag(— 1, 1, 1, 1) is the Minkowski met- 
ric and la is an outgoing null vector with respect to both 
the Minkowski and the Schwarzschild metric; the outgo- 
ing null vector written explicitly in terms of the coordi- 
nates (i, x, y, z) has the form 



Ir.^ \l 



x-xq y-yo z- zq 



(18) 



a 



1 + 2M/r ' 



and 



Pi = -^ (x -xo,y- yo, z - zq) . 



(19) 



(20) 



The computational parameters of the simulation are 
the same as in the previous section: N^ points with 
N — 361; outer boundary locations at [—2.5,-1-2.5] in 
the X, y, z directions; a Courant-Friedrichs-Lewy factor of 
A = 1/4; and artificial viscosity parameter of e = ft,^/16. 
However, Ki 550 evolution steps to the t = —1.910 sec 
time slice are needed in this particular case because of 
the shear of the null cone and the movements of the pul- 
sars, both of which will be discussed later in this section. 

As noted in Section |IV[ the vertex of the null cone is 
difficult to handle because of the finite resolution of the 
grid used to evolve the null cone into the past. In the 



fiat Minkowski space treated in Section IV the null cone 
can be described analytically, and is a shear-free met- 
ric sphere expanding from the vertex, with radius equal 
to the elapsed Minkowski time. Thus in Section IV we 
set the data at a Minkowski time that is 0.1 sec earlier 
(Ai = —0.1 sec) than the event we are coordinatizing — 
when the backward light cone sphere had a radius of 
0.1 sec. This sphere is sufficiently large so that our dis- 
cretization adequately resolves it, and as we follow the 
expanding light cone into the past, the relative resolu- 
tion becomes even better. In the Minkowski case this 
data setting method contributes to minimal total errors 
(a few times the discretization size). 

We use a similar method to initialize the null "cone" 
here in the black hole spacetime. (The null cone here is 
not a spherical cone because of shear due to the pres- 
ence of the black hole.) Again we analytically set the 
boundary condition at a small coordinate time (0.1 sec) 
in the past of the event we are coordinatizing. To pro- 
vide a technique for general spacetime, we define an ap- 
proximate method that will work in any spacetime: we 
evaluate the metric at the event being coordinatized, and 
assume that it is constant in the small region needed to 
propagate the null surface backward for a small arbitrary 
coordinate time. 

In the current example, we will find the emission co- 
ordinates of the event given by {t,x,y,z) = (0.1,0,0,0) 
in the Eddington-Finkelstein coordinates where the black 
hole is offset to {xo,yo,ZQ) = (2.5,0,0) and has a mass 
M = 0.25. Therefore, 2M/r = 1/5 at (0.1,0,0,0), 
and from Eqs. ^jl, (18) the metric at that point is 
9u = -4/5; gtx = -1/5; Oxx = 6/5; gyy = g^z = 1- 

One proceeds by choosing a small interval of 
Eddington-Finkelstein coordinate time At (here —0.1 
sec) and then solving for the values of the 3-space points 
{Ax, Ay, Az) that satisfy 



In Eddington-Finkelstein coordinates, the quantities 



= guAf + 2gtxAtAx + gijAx'Ax^ 



(21) 



using the constant values of the metric coefficients de- 
fined above. The 3-space points found define the shape 
of the t = constant — At shce of the backward nuh cone, 
which initializes the eikonal data for further evolution 
into the past. 

Contrary to the Minkowski case, now the pulsars 
are freely falling toward the black hole located at 
{xo,yo,zo) — (2.5,0,0). As an initial condition we de- 
mand that at time t = sec, all five pulsars are located 
at coordinate positions with the same values as used in 
the previous section, with zero velocity. However, due 
to black hole gravitational acceleration, the pulsars are 
moving at other times. Because they begin at rest in the 
black hole frame, the pulsars have only a radial velocity — 
and no angular velocity — toward the black hole. There- 
fore, to measure the emission coordinates of an event at 
the origin, we need to first understand the radial geodesic 
motion of the freely-falling pulsars. 

The Schwarzschild metric written in spherical 
Eddington-Finkelstcin coordinates is 



ds^ = - f 1 - ^^ ) df2 



UI 



2M\ 
r J r 



drdt 



1 



^^d.^ 



r I 



(22) 



and the Schwarzschild metric written in the 
Schwarzschild coordinates is 



ds' 



2M\ ,-2 
d? 



r J 



^V'd.^ 



r / 



-f r^(dr +sin-'6'd(y9^). 



(23) 



The bar on t here is used to distinguish the Schwarzschild 
coordinate time t from the Eddington-Finkelstein coor- 
dinate time t. 

From the expressions of Schwarzschild metric in the 
two coordinate systems, it is clear that the area co- 
ordinate r in Eddington-Finkelstein coordinates is the 
same area coordinate r in Schwarzschild coordinates. 
This guarantees that a pulsar's radial geodesic motion 
found using Schwarzschild coordinates, when expressed 
in terms of the pulsar proper time r, will have identi- 
cally the same expression in Eddington-Finkelstein coor- 
dinates. Let us proceed in the Schwarzschild coordinates. 

Purely radial geodesic motion of a freely falling pulsar 
can be described by the parametric equations ^22j: 



R 



(1 -I- cos 77), 



(24) 



2 yiM) 



('7 + sinr/), 



(25) 



and 



2M 



V2Af 



1/2- 



R ( R 

2" V2M 



1/2 



sm?7 



2M In 



(i?/2Af- 1)1/2+ tan (r;/2) 
(i?/2M- 1)1/2- tan (r;/2) 



(26) 



r 



where r\ is the parameterization, and R is the apastron — 
the areal distance at which the pulsar has zero velocity. 
Note that both t = and r = when ry = 0. As a specific 
choice in our simulation, we also specify that r = 0, for 
each pulsar, occurs at Eddington-Finkelstein coordinate 
time i = 0. 

We have obtained the expression for Schwarzschild 
time, but the Schwarzschild metric we supplied to the nu- 
merical simulation is expressed in Eddington-Finkelstein 
coordinates. Thus, we need to relate the Eddington- 
Finkelstein time and the Schwarzschild time, which can 
be done by relating the spacetime interval in Eqs 
and ( 23 1 to obtain 



(22) 



dp 



At 



2M/r 



Ar 



(27) 



1 - 2M/r 
Taking the positive root of this equation and integrating. 



we obtain the expected expression: 

t = t + 2Afln(r-2M) + C 



(28) 



where C is an arbitrary constant that depends on the 
initial condition. In our simulation, the initial condition 
is at t = t = 0, r = i?. Substituting this and solv- 
ing for C, we obtain the final equation that relates the 
Schwarzschild time i with the Eddington-Finkelstein time 



t = t + 2M\n{r- 2M) - 2M In {R - 2M). 



(29) 



Notice that R depends on which pulsar is under consider- 
ation, so we are setting data that have t = Q — constant, 
but have a different t for each pulsar. We now have the 
complete relations needed to describe the purely radial 
geodesic motion of the pulsars. 

As we run our simulation to obtain the null surface 
r, the pulsars move along the geodesies as described by 



Eqs. ([24]), ([26]), and (|29|. Figure [2] depicts the intersec- 
tions of the past light cone generated by the simulation 
at the equator, z — 0. The world lines of the first four 
pulsars are again depicted by four bold lines and their 
intersections with the null cone by dots. 

The curvature of pulsar #2's world line is most evident 
in Figure [2] Note also in Figure [2] the shape of the light 
cone: the light cone moves further away from the black 
hole as we go further into the past, as expected from 
intuition. (It is falling into the black hole as time goes 
forward.) The cone also flattens out as we go further 
into the past due to the tidal field of the black hole. 
These movements of the light cone are a main reason 
why the measured emission coordinates in Table [iT] the 
results of the measurements in Schwarzschild geometry, 
differ from the results in Table [l] The proper times r" 
in Table [ll] could in principle be obtained by integration 
back along the pulsar world lines. However, with the 
analytic results, r" are in fact obtained by computing 77 
for each pulsar at the event at which it intersects the null 
cone from Eqs. (24), (26), and (291, and then obtaining 
T from Eq. ([25]). 



Although we do not develop a complete analytical so- 
lution for the photons in this case (see [11] for such a 
solution), the configuration is qualitatively comparable 
to the Minkowski case discussed earlier. Therefore, we 
assume a similar level of error in the Schwarzschild case 
and indicate that in Table HTl 



Pulsar 


# 


r 


t 


X 


y 


z 


1 




-0.364 ±0.003 


-0.399 ±0.003 


-0.498 


0.000 


0.000 


2 




-1.432 ±0.003 


-1.840 ±0.003 


1.170 


0.000 


0.000 


3 




-0.641 ± 0.003 


-0.715 ±0.003 


0.007 


-0.748 


0.000 


4 




-1.112 ±0.003 


-1.233 ±0.003 


0.017 


1.241 


0.000 


5 




-0.618 ±0.003 


-0.701 ±0.003 


0.309 


0.398 


0.498 



TABLE II. Results of emission coordinates of an event X = 
{t, X, y, z) — (0.1, 0, 0, 0) in Eddington-Finkelstein coordinates 
in Schwarzschild spacetime. The pulsars are freely falling to- 
ward a static black hole located at {xo,yo,zo) — (2.5,0,0). 
t denotes the Eddington-Finkelstein times when the pulsar 
world lines intersect with the observer's past light cone; at 
i = 0, all pulsars have zero velocity, r lists the proper times: 
the emission coordinates of the reception event K can be 
recorded as the collection of any four of the t°'. x, y, and z are 
the Eddington-Finkelstein spatial coordinates of the pulsars 
when the intersections occur. 



necessary for the stability of this curved space simulation. 
The eikonal evolution suffered abrupt failures after ~ 200 
integration steps; reinitializing every 75 integration steps 
controlled this behavior. Results were independent of 
the frequency of reinitialization, so long as it was no less 
frequent than every 75 steps. 



VI. ASTROPHYSICAL APPLICATION: FIXING 
GAUGE 




FIG. 2. Rays: the past null cone F of the measurement 
event 3?. We use Eddington-Finkelstein coordinates in a 
Schwarzschild spacetime. Star marks the location of event 
3? = {t, X, y, z) = (0.1, 0, 0, 0) that we are coordinatizing. Bold 
lines: world lines of four z = pulsars freely falling toward a 
stationary black hole located to the right of the volume plot- 
ted at {xo,yo,ZQ) = (2.5,0,0). Dots mark the intersections 
between the null surface and the pulsar world lines. The co- 
ordinates of these intersections are recorded in Table [III See 
the text for a discussion of the gap between K and the top of 
the (truncated) null cone. 



The numerical reinitialization described in Eq.( 16 1 was 



Any set of four proper emission times may be used as 
emission coordinates, which we generically call C". Since 
here we collect pulses from five pulsars, there are five such 
sets and we introduce a numbering ^C" over all coordinate 
systems created by combining selected t^ . A typical set 
in our case might be: 2C" = 1''"^: ■'"'^:'''^,t^}; we restrict 
the possible choices by requiring the numbering to be 
ordered, increasing, in the quadruple. We note again 
that |10l shows how to combine results from more than 
four pulsar sources in a way that produces coordinates 
with reduced uncertainties. 

Since we are dealing with a particular event D? in a 
given spacetime, all the iC" ^re coordinate transforma- 
tions of one another. Since the source pulsars are inde- 
pendent, the transformations are discontinuous: no con- 
tinuous path of small transformations joins them. But 
besides these transformations, once a set of source pul- 
sars is decided on, there is another group of continuous 
gauge transformations, which we will discuss below. 

The emission coordinates — say {t^ , t"^ , t^ , t'^) — form 
a quadruple of proper emission times of the measured 
pulses. For now we assume that the pulses are closely 
enough spaced such that each t" can be viewed as a 
continuous time signal determined to arbitrary preci- 
sion. (Our simulations made this assumption and we 
used other means — not pulse counting — to determine the 
source proper times.) In practice one would interpolate; 
the interpolation process will be discussed briefly below. 



It is clear that this coordinate system has a gauge 
group: afSne transformations on each of the r" — >■ 
A"r" + B" (no sum on a), where A" and _B" are fi- 
nite and we restrict A" to be positive. Clearly B" is an 
offset (e.g. Eastern time vs Pacific time) and A" is a 
clock rate factor, or could be viewed as a function of the 
time unit chosen (e.g. seconds vs hours). 

To construct a consistent coordinate system, A" and 
B" must be chosen and held fixed for all reception events. 
(In our simulations we set A" = 1, and defined B" by 
demanding that t" — when the pulsar coordinate time 
t = 0.) 

Setting the gauge is intertwined with other steps in 
establishing the emission coordinates. For practical pur- 
poses, one may choose to construct a consistent coordi- 
nate system, using one or more central master stations, 
by proceeding as follows: 

1. Accumulate data on potential source pulsars. Com- 
pare pulse arrival stability against the best atomic 
clock. This will require the removal of known detec- 
tor motion in the Solar System, and further poly- 
nomial fitting of the pulse arrival times. Transfer 
time of arrival solutions to a fiducial point, such as 
the barycenter of the Solar System. 

2. Using good atomic clocks which hold stability over 
many pulse periods, interpolate the intervals be- 
tween source pulses. 

3. Define t*^ = to correspond to the arrival of a 
specific pulse from each source (3 at the receiver. It 
is intended that this is done at some finite specific 
time, while the pulse stability is being observed. 

Steps 1 and 2 provide precise interpolation into the 
intervals between source pulses, and partly imply a gauge 
(defining A") for the coordinate system. Step 3 defines 
the origin of the emission coordinate system, including 
fixing the offset _B". As the receiving station 3? ages and 
moves, its emission coordinates will move in a smooth 
way. 



being implemented to improve our code accuracy. The 
current code employs second-order discretization. This 
is being improved stepwise to fourth and ultimately to 
eighth order discretization. The current code is a uni- 
grid code, so the same numerical error limit applies at all 
points of the grid. The backward wavefront (backward 
null cone) curvature is greatest near the event 31, but is 
small for most of the evolution from sources to 31. We will 
implement multiresolution so that the regions where the 
null cone is most highly curved are well resolved. Even 
with this improvement, we still require setting data on 
a small "sphere" just to the past of the event 31. Care 
will be taken that this data setting is consistently con- 
vergent with the rest of the computation. Additionally 
we are studying the behavior of the eikonal differential 
equation, to optimize for accuracy its translation into 
computational terms. 



VIII. CONCLUSION 

We have presented a robust numerical method to mea- 
sure the emission coordinates of an event in any generic 
spacetime configuration. Our method uses a computa- 
tional evolution of the eikonal equation describing the 
backward light cone. 

We applied our method in two different numerical 
simulations: one in Minkowski spacetime and one in 
Schwarzschild spacetime. In both simulations, we found 
that our methods are reliable in measuring the emis- 
sion coordinates. Errors in the measurements can be 
attributed to the numerical accuracy of the simulations, 
mainly due to the interpolation routine and the resolu- 
tion of the three dimensional finite difference mesh. We 
anticipate that with a higher resolution, we will be able 
to reduce the errors in our calculation of the emission 
coordinates. 

Although these numerical simulations are preliminary, 
the same method can be used for pertinent problems, 
such as measuring the emission coordinates of the Earth, 
and therefore the Earth's trajectory. 



VII. COMPUTATIONAL ACCURACY 
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